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Abstract:  Effects  of  uniform  temperature  and  linear  temperature  gradi¬ 
ent  on  flow  rate  and  density  profile  of  fluid  driven  from  a  source  to  an  open 
system  is  studied  by  a  Monte  Carlo  simulation  in  three  dimensions.  The 
steady-state  density  profile  with  uniform  temperature  differs  significantly 
from  that  of  the  temperature  gradient  at  all  but  high  temperature  regimes 
where  the  profile  is  a  linear  density  gradient;  the  crossover  to  linear  density 
gradient  of  the  profile  is  sensitive  to  range  of  temperature  variation.  The 
response  of  the  flow  rate  density  to  the  temperature  is  nonlinear. 
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Understanding  the  flow  of  methane  gas  below  the  ocean  floor  is  of  consid¬ 
erable  interest  particularly  in  context  to  geochemical  and  geophysical  charac¬ 
teristics  of  methane  hydrate  in  marine  sediments  [1-10].  At  the  temperatures 
and  pressures  within  the  marine  sediments,  the  gas  (usually  Methane)  is  in 
solution.  Systematic  experimental  observations  of  flow  and  density  profile 
of  methane  gas  are  severely  limited  due  to  uncontrollable  parameters  such 
as  changes  in  pressure  and  temperature  gradients  and  morphological  varia¬ 
tions  of  the  host  matrix,  i.e.,  the  evolution  of  porous  ground  [5-7].  One  of 
the  major  problems  in  modeling  such  a  complex  system  is  to  incorporate  all 
the  details  on  a  large  scale  [11].  Therefore,  identifying  their  contributions, 
such  as  effects  of  temperature  gradient,  pressure  gradient,  porosity  etc.  on 
such  properties  as  flow  rate  of  methane  and  its  density  profile  is  very  difficult 
[12,13],  Obviously,  this  requires  long  term  study  involving  model  develop¬ 
ments  and  analysis  of  resulting  data.  It  would  be,  therefore,  desirable  to 
focus  on  one  parameter,  i.e.,  temperature  gradient,  at  a  time  and  under¬ 
stand  its  role,  as  attempted  in  this  paper  before  including  other  parameters 
in  the  model. 

We  consider  a  discrete  lattice  of  size  Lx  Lx  Lz.  As  in  our  earlier  model, 
the  bottom  plane  (z  =  0)  is  connected  to  a  source  of  fluid  (i.e.,  methane 
in  solution)  modeled  by  an  interacting  lattice  gas.  The  top  plane  (z  =  Lz) 
of  the  lattice  is  open,  i.e.,  a  gas  particle  can  escape  from  top.  In  order  to 
consider  interaction,  we  introduce  an  interaction  density  (pi)  to  each  lattice 
site  i  :  pi  =  1  if  site  i  is  occupied  by  a  fluid  particle  and  pi  —  —1  if  site  i  is 
empty.  The  interaction  Hamiltonian  is  given  by, 

E  =  U^2  PiPj,  (1) 

ij 

where  U  is  the  interaction  strength;  we  set  U  =  1  in  unit  of  Boltzmann 
constant  ( kB )  but  U  can  be  varied  if  needed  in  order  to  incorporate  specific 
effects  such  as  drag/viscosity  of  the  host  matrix.  The  summation  is  eq.  1 
is  restricted  to  nearest  neighbor  sites.  Initially,  the  entire  lattice  is  empty 
except  the  first  plane  (z  —  1)  connected  to  the  source  which  is  filled  with 
fluid,  one  particle  at  each  site. 

Two  thermal  systems  are  considered:  (i)  homogeneous  temperature  (T) 
throughout  the  sample  from  bottom  to  top  and  ( ii )  a  linear  temperature 
gradient, 

T(z)  =  (Tb  -  Tt)/z,  (2) 
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where  T(z )  is  the  temperature  of  2— th  plane  and  Tb  =  T(z  —  1),  T)  =  T(z  = 
Lz).  The  temperature  at  the  top  plane  T)  =  0  is  fixed  while  the  temperature 
at  the  bottom  plane  Tb  is  varied.  Fluid  particle  executes  its  stochastic  motion 
with  the  Metropolis  algorithm:  a  particle  at  site  i  and  one  of  its  neighboring 
site  j  are  chosen  randomly.  If  site  j  is  empty  then  we  accept  the  move  with  the 
Boltzmann  probability  exp(-AE/kBT),  where  A E  is  the  change  in  energy 
due  to  corresponding  configurational  changes,  i.e.,  from  pi  =  1,  Pj  =  —1  to 

=  -l,pj  =  1.  Energy  and  temperature  are  measured  in  arbitrary  unit 
(normalized  by  the  Boltzmann  constant  ( kB ))•  Periodic  boundary  condition 
is  used  along  x-  and  y-directions  and  open  condition  along  (z-)direction. 
Thus,  a  fluid  particle  can  escape  from  the  top  (z  =  Lz)  but  cannot  go  below 
the  bottom  layer  {z  =  1).  An  attempt  to  move  each  particle  once  is  defined  as 
one  Monte  Carlo  Step  (MCS).  As  soon  as  a  site  in  the  bottom  plane  (z  =  1) 
is  vacated  by  a  mobile  gas  particle,  it  is  filled  by  a  particle  from  the  source. 
Note  that  the  number  of  particles  is  not  conserved. 

As  the  simulation  proceeds,  the  fluid  flows  from  bottom  to  top  and  the 
density  profile  evolves.  During  the  simulation  we  monitor  the  overall  mobility 
of  fluid  particles  and  their  center  of  mass,  flow  rate  of  fluid,  and  their  density 
profile.  The  simulation  is  repeated  for  a  large  number  of  independent  runs 
in  order  to  estimate  the  average  values  of  the  steady-state  flow  rate  and 
equilibrium  density  profile.  In  the  steady-state,  the  flow  rate  of  fluid  is 
independent  of  time  and  the  density  profile  becomes  constant  (see  below). 
The  simulation  is  carried  out  on  different  lattice  sizes  but  most  simulations 
are  performed  with  303  and  603  samples;  we  have  not  seen  severe  finite  size 
effects  on  the  qualitative  behavior.  Figure  1  shows  the  variation  of  net  mass 
of  the  fluid  flow,  i.e.,  the  fluid  flux  (4>)  with  MCS  time  along  the  vertical  (z) 
direction.  The  linear  derivative  of  mass  flow  with  time  leads  to  fluid  current, 

*=*•  (3) 

and  current  density, 

j  =  i/A,  (4) 

where  A  =  L2,  the  cross-sectional  area.  The  linear  variation  of  the  data 
suggests  that  the  flow  rate  is  independent  of  time.  This  implies  that  our 
system  has  reached  the  steady-state. 

The  response  of  the  flow  rate  density  to  temperature  is  nonlinear  (figure 
2).  Generally  the  flow  rate  response  is  low  at  low  temperatures  and  high 
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temperatures  with  a  large  change  in  response  rate  in  between.  The  nonlinear 
pattern  of  the  response  with  uniform  temperature  is  not  much  different  from 
the  linear  temperature  gradient  except  the  range  of  temperatures. 

The  steady-state  density  profiles  at  different  temperatures  (case  ( i ))  and 
temperature  gradients  (case  (ii))  are  presented  in  figure  3.  We  note  a  signif¬ 
icant  difference  in  the  density  profiles  in  systems  with  uniform  temperature 
and  linear  temperature  gradients.  With  the  uniform  temperature,  a  sig¬ 
moidal  form  of  profile  develops  except  at  high  temperatures  (T  =  2.0, 3.0) 
where  the  density  decays  approximately  linearly  from  bottom  to  top.  At  low 
to  moderate  temperatures  ( T  =  0.0167  -  1.0),  three  linear  density  gradients 
are  observed:  a  slow  linear  decay  from  bottom  in  lower  half  of  the  sample 
[z  =  1  -  12)  is  followed  by  a  large  density  decay  in  few  planes  (z  =  15  -  19) 
before  the  onset  of  a  low  linear  decay  toward  the  top.  In  presence  of  lin¬ 
ear  temperature  gradient  on  the  other  hand,  the  approach  to  linear  density 
profiles  at  high  values  of  Tb  (=3.0  -  10.0),  is  more  step-like  (see  fig  3(b)); 
the  crossover  to  linear  decay  of  density  in  upper  part  of  the  sample  is  rather 
slow.  In  either  case,  the  effect  of  temperature  change  is  very  important  in 
orchestrating  the  shape  of  the  density  profile.  Note  that,  in  addition  to  tem¬ 
perature,  the  density  gradient  also  plays  a  significant  part  in  driving  the  fluid 
density  from  bottom  to  top.  Thus,  if  the  measurement  of  temperature  profile 
below  the  ocean  floor  is  available,  i.e.,  uniform  and/or  linear  gradient,  then, 
one  would  be  able  to  gain  an  insight  into  the  density  profile  of  methane  in 
solution  and  possibly  the  location  of  methane  hydrate  formation. 

Connecting  such  qualitative  results  to  quantitative  field  measurement  is 
one  of  the  major  problem  in  such  simulations.  Scaling  and  calibration  of  the 
physical  quantities  and  variables  (time  steps,  temperature,  and  lattice  size) 
with  the  measured  quantities  are  needed  in  order  to  gain  insight  into  the 
applicability  of  these  findings  in  the  range  of  variables.  Unfortunately,  field 
measurements  for  systematic  variation  of  the  density  profile  or  the  flux  rate  of 
methane  [13]  with  depth  and  temperature  gradient  is  not  available  at  present. 
Neither  it  is  feasible  to  approach  the  relevant  scale  of  marine  geosystem  in 
computational  models  from  microscopic  scale.  Therefore,  it  is  rather  difficult 
to  connect  these  results  to  limited  field  measurements  in  a  meaningful  way. 
However,  it  is  worth  pointing  out  the  order  of  scales  even  with  our  coarse¬ 
grained  approach.  Typically,  a  temperature  gradient  below  the  ocean  floor 
is  about  20°C/km.  Our  temperature  gradient  A T/z  ~  20°C/km,  leads  to 
energy  of  gas  of  the  order  of  10-19  J  ~  1  ev  if  we  scale  the  temperature  T  with 
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the  Boltzmann  constant  and  z  ~  1km.  Equipartition  of  energy,  provides  an 
estimate  of  particle  speed  ~  37  m/s.  The  variation  in  the  the  temperature 
gradients  are  selected  to  explore  the  changes  in  density  profile  and  response 
of  the  flow  rate.  Now,  let  us  look  at  the  advance  of  the  fluid  front  to  estimate 
the  time  steps.  A  typical  variation  of  the  center  of  mass  of  the  fluid  with 
time  is  presented  in  figure  4.  This  gives  the  front  velocity  as  0.02  lattice 
constant  /  time  steps.  Since  we  have  assume  the  size  of  the  sample  to  be  1.0 
cm,  one  lattice  constant  is  0.033  cm  on  a  303  lattice.  In  order  to  match  the 
fluid  flow  rate  of  1  mm/yr  [13]  in  field  measurement,  our  time  step  should  be 
of  the  order  of  a  month.  However,  if  the  field  measurement  for  the  flow  rate 
is  lower  then  our  MC  step  has  to  be  increased  accordingly.  Note  that  we  do 
not  have  field  measurements  suitable  for  our  idealized  model  and  therefore, 
these  estimates  are  mainly  to  illustrate  the  method. 

In  summary,  a  computer  simulation  model  is  presented  to  show  that  the 
density  profile  of  fluid  driven  from  a  source  into  an  open  system  depends  on 
temperature  and  the  temperature  gradient:  specific  profiles  are  proposed  for 
a  range  of  temperature  and  temperature  gradients.  The  flow  rate  density 
responds  nonlinearly  to  variations  in  temperature  and  gradient.  Even  the 
simple  ’’toy  model”  like  this  one  is  computationally  intensive  for  the  complex 
issues  addressed  here  but  provides  insight  into  evolution  of  density  profile 
and  flow  rate  that  would  be  encountered  with  methane  in  solution  moving 
through  marine  sediments. 
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Figure  Captions: 

Figure  1:  Net  mass  of  gas  flow  along  z-direction  versus  MCS  time.  603 
sample  is  used  with  10  independent  runs  at  TJ,  =  0.5 

Figure  2:  Variation  of  the  flow  rate  density  of  gas  with  temperature.  Sample 
of  size  303  is  used  with  up  to  50  independent  runs. 

Figure  3:  Density  profile  for  uniform  temperature  (T)  (a)  and  linear  tem¬ 
perature  gradient  (b)  with  various  T  and  T&;  samples  of  size  303  is  used  with 
20-50  independent  runs. 

Figure  4:  Variation  of  the  center  of  mass  of  the  gas  with  time  step  on  a  303 
lattice  for  three  temperature  gradients.  Statistics  is  the  same  as  in  figure  3. 
The  dotted  lines  are  the  least  square  fit  and  give  a  slope  of  about  0.02. 


7 


CCTT'-CO 

t-COOCO 


OO-COUTCOO 

dddddr-Ovbd 

ll  II  ii  II  ll  II  II  II 


CD  t3" 

o  d 


Ajisuen  0uein0|/\| 


